function value = subsref(self,index)
%*=+--+=#=+--      UAS Severe Weather Simulation Softeware        --+=#=+--+=#*%
%          Copyright (C) 2013 Regents of the University of Colorado.           %
%                             All Rights Reserved.                             %
%                                                                              %
%    This program is free software: you can redistribute it and/or modify      %
%    it under the terms of the GNU General Public License version 2 as         %
%    published by the Free Software Foundation.                                %
%                                                                              %
%    This program is distributed in the hope that it will be useful,           %
%    but WITHOUT ANY WARRANTY; without even the implied warranty of            %
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             %
%    GNU General Public License for more details.                              %
%                                                                              %
%    You should have received a copy of the GNU General Public License         %
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.     %
%                                                                              %
%           Jack Elston                                                        %
%           elstonj@colorado.edu                                               %
%                                                                              %
%*=+--+=#=+--                 --+=#=+--+=#=+--                    --+=#=+--+=#*%

  switch index.type
    case '()'
      time=index.subs(1);
      time=time{1};

      if(size(index.subs,2) == 2) 
        pos=index.subs(2);
        pos=pos{1};

        pos(1) = pos(1) - self.simdata.frame.u*(time-self.simdata.t0);
        pos(2) = pos(2) - self.simdata.frame.v*(time-self.simdata.t0);

      else 
        x=index.subs(2);
        x=x{1};
        y=index.subs(3);
        y=y{1};
        z=index.subs(4);
        z=z{1};

        x = x - self.simdata.frame.u*(time-self.simdata.t0);
        y = y - self.simdata.frame.v*(time-self.simdata.t0);

        x = floor((x - self.simdata.x0)/self.simdata.xs) + 1;
        y = floor((y - self.simdata.y0)/self.simdata.ys) + 1;
        %z = floor((z - self.simdata.z0)/self.simdata.zs) + 1
				z = min(find(self.simdata.z >= z)) % work around to be able use non-uniform grids
        t = floor((time - self.simdata.t0)/self.simdata.ts) + 1;
      
        value.u    = self.simdata.u(x,y,z,t) + self.simdata.frame.u;
        value.v    = self.simdata.v(x,y,z,t) + self.simdata.frame.v;
        value.w    = self.simdata.w(x,y,z,t);
        value.rain = self.simdata.rain(x,y,z,t);
        value.hail = self.simdata.hail(x,y,z,t);

        return
      end
      
      x = self.simdata.x;
      y = self.simdata.y;
      z = self.simdata.z;
      t = self.simdata.time;

			%FIXME -- need to validate values out with data
      if self.simdata.znum == 1 && self.simdata.tnum == 1
        u    = interpn(x,y,self.simdata.u,   pos(2),pos(1));
        v    = interpn(x,y,self.simdata.v,   pos(2),pos(1));
        w    = interpn(x,y,self.simdata.w,   pos(2),pos(1));
        rain = interpn(x,y,self.simdata.rain,pos(2),pos(1));
        hail = interpn(x,y,self.simdata.hail,pos(2),pos(1));
      elseif self.simdata.znum == 1
        value = 0;
        return %not handled yet
      elseif self.simdata.tnum == 1
        value = 0;
        return %not handled yet
      else
        %--- constraint --%
        c(1) = floor((pos(1) - self.simdata.x0)/self.simdata.xs) + 1;
        if (c(1) <= 1)  c(1) = 2; end;
        if (c(1) >= self.simdata.xnum)  c(1) = self.simdata.xnum -1; end;
        c(2) = floor((pos(2) - self.simdata.y0)/self.simdata.ys) + 1;
        if (c(2) <= 1)  c(2) = 2; end;
        if (c(2) >= self.simdata.ynum)  c(2) = self.simdata.ynum -1; end;
        %c(3) = floor((pos(3) - self.simdata.z0)/self.simdata.zs) + 1;
				c(3)  =min(find(self.simdata.z >= pos(3))); % work around to be able use non-uniform grids
        if (c(3) <= 1)  c(3) = 2; end;
        if (c(3) >= self.simdata.znum)  c(3) = self.simdata.znum -1; end;
        c(4) = floor((time - self.simdata.t0)/self.simdata.ts) + 1;
        if (c(4) <= 1)  c(4) = 2; end;
        if (c(4) >= self.simdata.tnum)  c(4) = self.simdata.tnum -1; end;

        x = x(c(1)-1:c(1)+1);
        y = y(c(2)-1:c(2)+1);
        z = z(c(3)-1:c(3)+1);
        t = t(c(4)-1:c(4)+1);
        
        subset.u = self.simdata.u(c(2)-1:c(2)+1,c(1)-1:c(1)+1,c(3)-1:c(3)+1,c(4)-1:c(4)+1);
        subset.v = self.simdata.v(c(2)-1:c(2)+1,c(1)-1:c(1)+1,c(3)-1:c(3)+1,c(4)-1:c(4)+1);
        subset.w = self.simdata.w(c(2)-1:c(2)+1,c(1)-1:c(1)+1,c(3)-1:c(3)+1,c(4)-1:c(4)+1);
        subset.rain = self.simdata.rain(c(2)-1:c(2)+1,c(1)-1:c(1)+1,c(3)-1:c(3)+1,c(4)-1:c(4)+1);
        subset.hail = self.simdata.hail(c(2)-1:c(2)+1,c(1)-1:c(1)+1,c(3)-1:c(3)+1,c(4)-1:c(4)+1);
				
        u    = interpn(x,y,z,t,subset.u,   pos(1),pos(2),pos(3),time);
        v    = interpn(x,y,z,t,subset.v,   pos(1),pos(2),pos(3),time);
        w    = interpn(x,y,z,t,subset.w,   pos(1),pos(2),pos(3),time);
        rain = interpn(x,y,z,t,subset.rain,pos(1),pos(2),pos(3),time);
        hail = interpn(x,y,z,t,subset.hail,pos(1),pos(2),pos(3),time);
        %--- constraint --%

%        u    = interpn(x,y,z,t,self.simdata.u,   pos(1),pos(2),pos(3),time);
%        v    = interpn(x,y,z,t,self.simdata.v,   pos(1),pos(2),pos(3),time);
%        w    = interpn(x,y,z,t,self.simdata.w,   pos(1),pos(2),pos(3),time);
%        rain = interpn(x,y,z,t,self.simdata.rain,pos(1),pos(2),pos(3),time);
%        hail = interpn(x,y,z,t,self.simdata.hail,pos(1),pos(2),pos(3),time);
      end

      u = u + self.simdata.frame.u;
      v = v + self.simdata.frame.v;
      
      value = [u,v,w,rain,hail];
    case '.'
      switch index.subs
        case 't0'
          value = self.simdata.t0;
        case 'ts'
          value = self.simdata.ts;
        case 'tf'
          value = self.simdata.tf;
        case 'tnum'
          value = self.simdata.tnum;
        case 'time'
          value = self.simdata.t0:self.simdata.ts:self.simdata.tf;
        case 'x0'
          value = self.simdata.x0;
        case 'xs'
          value = self.simdata.xs;
        case 'xf'
          value = self.simdata.xf;
        case 'xnum'
          value = self.simdata.xnum;
        case 'x'
          value = self.simdata.x;
        case 'y0'
          value = self.simdata.y0;
        case 'ys'
          value = self.simdata.ys;
        case 'yf'
          value = self.simdata.yf;
        case 'ynum'
          value = self.simdata.ynum;
        case 'y'
          value = self.simdata.y;
        case 'z0'
          value = self.simdata.z0;
        case 'zs'
          value = self.simdata.zs;
        case 'zf'
          value = self.simdata.zf;
        case 'znum'
          value = self.simdata.znum;
        case 'z'
          value = self.simdata.z;
        case 'frame'
          value = self.simdata.frame;
        case 'source'
          value = self.settings.dataset;
        case 'settings'
          value = self.settings;
        case 'simdata'
          value = self.simdata;
        case 'offset'
          value = self.simdata.offset;
      end
  end
